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Extended Abstract 


This paper will present a complete formulation of the optimal control problem for at- 
mospheric ascent of rocket powered launch vehicles subject to usual load constraints and 
final condition constraints. We shall demonstrate that the classical finite difference method 
for two-point-boundary- value-problems (TPBVP) is suited for solving the ascent trajectory 
optimization problem in real time, therefore closed-loop optimal endoatmospheric ascent 
guidance becomes feasible. Numerical simulations with a the vehicle data of a reusable 
launch vehicle will be provided. 

1. Ascent Guidance Problem Formulation 

The equations of motion of the RLV in an inertial coordinate system can be expressed as 


V = g(r) + A /m(t) + Tl b /m(t ) + N/m(t) (2) 

where r and V are inertial position and velocity vectors; g the gravitational acceleration; 
T the thrust magnitude; A and N are aerodynamic axial and normal forces, respectively; 
1„ the unit vector defining the RLV body longitudinal axis; m is the mass of the RLV. The 
magnitudes of the aerodynamic forces and thrust are modeled by 


A = \pV? 5 re /CU(Mach, a) 

A 

(3) 

N = ^pV?S T efC N ( Mach, a) 

(4) 

T = T vac + AT(r) 

(5) 


where p is the atmospheric density, and V r is the magnitude of the Earth-relative velocity 
V r = v — u>e x r with &e being the Earth angular rotation rate vector. The axial an 


‘Associate Professor, Department of Aerospace Engineering and Engineering Mechanics, Associate Fellow 

AIAA. Email: plu@iastate.edn _ . . , . 

t Graduate Research Assistant, Department of Aerospace Engineering and Engineering Mechanics 


1 


normal aerodynamic coefficients C A and C N are functions of angle of attack a and Mach 
number. They are usually expressed in analytical forms by curve-fitting tabulated data. 1 he 
reference area 5 re / is a constant. The vacuum thrust T vac may be time varying w en e 
thrust needs to be throttled down to meet axial acceleration constraint. The thrust loss 
inside the atmosphere AT < 0 is a function of altitude through the dependence of AT on 
ambient pressure. Mass m(t) is an explicit function of time. The mass flow rate wifi be 
reduced by the same percentage as the thrust when the thrust is throttled down. 

If the RLV symmetric plane is assumed to be always the plane formed by the body-axis 
1^ and the Earth-relative velocity vector V r , we can further have 

A = -Alb, N = Nl n ( 6 ) 


with the unit vector of the body 


normal axis l n defined by 


In 


(lftXlyJ 

11 ||1» X lv,ll 


where lvv = V r /V^. The angle of attack is then 

cos a = lblv r 


(7) 

( 8 ) 


The following expression for l n is valid for both o > 0 and a < 0 

l„ = l*xiiA^ W 

sma 

Note that in this formulation the sideslip angle (3 = 0 when there is no wind. 

The launch conditions for r 0 and V 0 are specified. The ascent guidance problem is to find 
the desired body-axis orientation l 6 (f) at each instant which determines the thrust direction 
and aerodynamic forces during atmospheric portion of the ascent. The final conditions wi 
be the engine-cutoff conditions which ensure insertion into the required orbit. These orbital 
insertion conditions can in general be written as k, 0 < k < 6, algebraic end conditions 

\P(r(t/), V(t/)) = 0, *eR k (10) 

In addition, there will be path constraints which limit the aerodynamic load and thrust 
acceleration of the RLV. These are expressed in terms of inequality constraints 

S(v,V,t)<0 (11) 

The mathematical tool used to solve this problem is the optimal control theory. In this 
setting a performance index is desired to be minimized, usually (but not necessarily ) the 
burn time of the rocket engine because it is directly related to the propellant usage. Denote 
the performance index by 

J = 


2 



The necessary conditions for the optimal control V b are given by defining the Hamiltonian 

H = p T V + p £ [g + (T — A)U/m + Nln/m] + A T S(r, V,t) + lb - 1) 

where /z is a scalar multiplier and A a vector multiplier of appropriate dimension. Then, the 
so-called costate equations and optimality condition for the optimal l* b are 


Pr = 


Pv = 


a H 

dr 

dH_ 

~dV 


^(Pr,Pv,r,V,lJ,t) = maxH(p r ,PK,r,V,l6,t) 


( 12 ) 

(13) 

(14) 


And the optimal solution satisfies the terminal constraints (10) and the transversality con- 
ditions 




d<S> 


(15) 

(16) 
(17) 


where v G R k is a constant multiplier vector. The first two conditions can be re^anged 
to yield 6 — k independent conditions involving only p / = (pj f Pv>) and x / = ( r / y ) * 
The general approach will be first finding the 6 - fc linear independent solutions of the 

homogeneous system ^ 

UJ*“° 

Let £(x/) € R 6 , i = 1,...,6 — fc be such solutions. Note that &’s are functions of X/. 
Transversality conditions (15) and (16) are then equivalent to 


(p/ + J^) T & = r i( p /> x /) = °> i = l,...,6-fc. 


(18) 


For a given problem, these conditions in above equation can often times be obtained more 
conveniently by using the terminal constraints (10) and taking dot products of Eqs. (15) 
and (16) with appropriate vectors related to the final state x/. 

When none of the inequality constraints in (11) are active, the multiplier vector A in the 
expression of H is zero. In such a case the condition (14) is equivalent to 


dH 

dl b 


= 0 


After much differentiation and algebraic operation, this condition can be shown to lead to 

lfc = CiPv + C2 V r ( 19 ) 
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where «’s are scalar functions of the state and costate. Therefore we conclude that the opti- 
mal body-axis lies in the plane formed by the primer vector Pv and relative velocity vector 
V . Note that as the atmospheric density decreases (approaching vacuum flight), c 2 -^U 
and Ci > 0. The optimal thrust vector becomes aligned solely with the primer vector Pv- 

The condition (19) suggests that the search for the optimal body-axis orientation can 
be reduced to a one-dimensional search in the plane of p v and V r . Let $ be the angle 
between the vectors p v and V r . Then it is straightforward to see that l b l Pv - cos(4> a) 
«nd 1 T 1 = sin(<& - a). Using these two relations in the expression of if, it is clear tnat 

maximizing H with respect to l b as in Eq. (14) is equivalent to dH/da = 0, which m turn 

rGSUltS m tan .($ - a) (T - A + N a ) -{A a + N) = 0 (20) 

with N a = dN/da and A a = dAjda. The above equation needs to be solved numerically 
for a (note that A, N , A a and N a are still functions of a). Once a is found, we have 

( sin a \ „ . r co s a - cos $ cos($ - a) 

sin 2 $ 


i; = 


Vsin^ 


1 PV + 




( 21 ) 


with l pv = p v/Pv- 

It should be stressed that for on-board guidance applications, it is imperative that great 
care be taken to determine the sign of $ in Eqs. (20) and (21). Not only the value and sign 
of a depend on the sign of $, but more importantly is the physical implication of the sign 
change of 4 > . Since the body y-axis can be defined by 

-i _ 1 Vr x jgv (22) 

y sin<fr 

When the vector l Vr x l pv changes direction by 180 degrees, the sign of $ should change. 
If the correct sign of 4> is identified when it becomes negative, the t/-body axis will have an 
instantaneous change of direction by 180 degrees. This would require an undesirable (may 
be even unfeasible) 180-degree roll of the vehicle. 

The typical path constraints (11) in an ascent guidance problem include aq, a, q, and 
acceleration limits. It can be shown that the condition (19) with different Ci and c 2 still 
holds true when these constraints become active. In these cases a is determined by the 
constraints rather than Eq. (20). The corresponding body axis l b is still found according to 
Eq (21) Hence the necessary conditions for the optimal control problem constitute a two- 
point-boundary-value problem (TPBVP) involving the state and costate equations (1-2) and 
(12-13), given initial launch conditions r 0 and V 0 , and the final conditions (10) and (18). 


2. Finite Difference Approach for Atmospheric Ascent Guidance 
2.1 Methodology 

Finite difference is one of the several techniques commonly used for TPBVPs. It tends 
to be more robust and less sensitive with respect to initial guesses as compared to shooting 
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methods. Suppose that the TPBVP at hand from the necessary conditions of the optimal 
control problem in the preceding section is put in the form of 

y = f{t, y) ( 23 ) 

#o(yo) = B fM = 0 ^ 

where y = (x T p T ) T € R 2n with n = 6. Let t/ be the specified final time. The 2 n 
boundary conditions in Eq. (24) axe from the given launch conditions, terminal constraints 
and transversality conditions. Note that the control l b has been expressed as a function of 
the y in the system equations through the optimality condition (21). To find the solution 
to the TPBVP, divide the time interval tf — to in to N subinterval of the same length 
h=(t f - t 0 ) /M. Let y fc = y(f 0 + kh) be the value of the solution at the node t k = t 0 + kh , 

fc = 0, . . . , M. At the middle point between t k - 1 and t k , denoted by t*- 1/2 = h — h/2, the 

differential equations (23) are approximated by central finite difference: 

^(yk - yik-i) = /( f k- 1/2) — 2 * 

Or equivalently, 

E k {y k ,y k -i) = yfc - Yk -1 ~ hf(h- 1 / 2 , y * y- ) = °> k = l,...,M -l. (26) 

In addition, the boundary conditions are denoted by 

£o(yo) = 5 0 (y 0 ) = 0 (27) 

= B f ( y/) = 0 (28) 

TVeat Y = (y^ yj ■ . • yJ{) T £ R 2n ( M+1 ) as the unknowns. The equal number of conditions 

316 E(Y) = 0 (29) 

where E = (E 0 Ex... E M ). Now the problem becomes a root-finding problem for a sys- 
tem of nonlinear algebraic equations. It has been rigorously established that under certain 
conditions on smoothness and the boundary conditions, the following holds true 4 

1. The original TPBVP and the finite difference problem have unique solution; 

2. The solution of the above finite difference problem y k is a second-order approximation 
to the solution of the TPBVP y *(t) at each t k , i.e., 

||y*(tfc) -yfc II = 0(h?) 


where lim^o 0(h?)/h 2 < 00 . 

For ascent guidance applications, since the time-to-go tf - 1 0 is decreasing, the accuracy 
of the finite difference solution will be higher and higher as h becomes smaller even for 
moderate number of nodes. 

2.2 Algorithm 
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The modified Newton method is probably the best suited algorithm for solving the prob- 
lem (29). Starting from an initial guess Y 0 , the search direction d, in the j-th iteration is 
determined by solving the linear algebraic equations 

d* = — E00-0, j = (30) 


dE(Yj-i) 

dY 


Then the update is given by 

Yj = Yj-i + crjdj, 0 < < 1, 


(31) 


The step size parameter (Tj is determined by the following criterion 

c Tj = maxj^-| E T (Yj_i + ^Tdj)E(Yj*_x + ^rdj) < E r (^i)E(l^„i)| (32) 


In other words, starting from Oj = 1, &j is halved repeatedly if necessary till the above con- 
dition is satisfied. 5 Convergence is achieved when ||E(Y i )|| < € for some preselected small 
£ > 0. The possible additional function evaluations required in checking the above step 
size condition pose negligible computational burden because function evaluations are not 
expensive in this setting. The result on the other hand is a much more robust algorithm, 
especially when the initial guesses are not close to the final solution. 


At the first glance, solving the linear system (30) may seem to be a formidable task 
because the dimension of the problem (2n(M + 1)) which can easily be over 1000. However, 
a close inspection reveals that the Jacobian matrix in (30) has a special sparse pattern due 
to the dependence of E^ on only y* and y*_i (cf. Eq. (26)). Therefore a very efficient 
algorithm, both in speed and storage, based on Gauss elimination and back substitution can 
be devised to solve the system (30). More details in implementation of such an algorithm 
are well documented in. 6 


The evaluation of the Jacobian dE/dY can certainly be done analytically. But we believe 
that simple forward finite differencing is more advantageous in this case. This is because 
unlike in a case where integrations of the differential equations are involved for each function 
evaluation, the function evaluation is fast. Using analytical Jacobian offers no computational 
speed benefits. On the other hand, analytical Jacobian will make the computer code signif- 
icantly more complicated because second-order partial derivatives of the right hand sides of 
the state equations are needed. Also, when some of the path constraints in Eq. (11) become 
active, the multiplier A will become functions of state and costate, adding more complexity 
to the analytical Jacobian. 

When using the finite difference algorithm, initial guess for the state as well as the 
costate are required at tjt, k = 0,1 The initial guess for x(t) can be obtained by 

linearly interpolate the given x 0 and an estimated x(f/) which is relatively easy to make 
because of the physical meaning of x(£/). For the costate p(£), & reasonable initial guess can 
be obtained from the costate history obtained by solving the problem with zero atmospheric 
density (vacuum flight). For vacuum flight, the finite difference algorithm converges in many 
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cases with even constant guess for p(t). Another possibility is to use the costate and state 
solutions from the analytical vacuum guidance algorithm which is described in the next 
section. The atmospheric solution is then obtained by applying homotopic continuation on 
the atmospheric density with 

p = rjp, 0 < r) < 1 

For each r?, p is used in place of atmospheric density in the aerodynamic forces terms. The 
homotopic parameter 77 starts at 0 (vacuum solution), and gradually increase to unity for 
full atmospheric solution. 

3. Analytical Vacuum Ascent Guidance 

While the finite difference approach in the preceding sections applies equally well to 
vacuum ascent guidance, a simpler and even better analytical approach exists. This approach 
combines a number of elegant results in optimal vacuum trajectory studies over the past three 
decades, and is summarized in Ref. 7 The key ingredients are the linear gravity simplification 
and the closed-form solution of the costate equation and nearly closed-form solution of the 
sate equation. At the beginning of each guidance update cycle, let r 0 be the position vector. 
The gravity acceleration g is then approximated by 

g = -^fr = -u, 2 r (34) 

r o 

where Pe is the gravitational parameter of the Earth. In an ascent guidance problem, the 
correct direction of the gravity is more important than the accuracy of its magnitude. This 
approximation preserves the change of direction of the gravitational acceleration with r. The 
magnitude of g will be slightly different from that of a Newtonian central gravity field. But 
when r 0 is continuously updated by the radius at beginning of each guidance cycle, the effect 
of this difference will be negligible. 


Let g 0 = pe/rl be the magnitude of the gravity acceleration a t r 0 . We normalize the 
equations of motion (1) and (2) with unit distance r 0 , unit time VWtfo, and unit velocity 
y/r^. The dimensionless equations of motion, with A = N = 0 for vacuum flight, become 

r = V (35) 

V = — r + T(r)l fc ( 36 ) 

where the Schuler frequency u ; has become unity in the normalized time, and T(t) = 
T va c/ m ( T )9o with T as the normalized time. Note that this normalization is done in each 
guidance cycle with the r 0 being the radius at the beginning of that cycle. The Hamiltonian 

now is ^ /o~\ 

H = p?V + p£[-r + T(r)l 6 ] + p{tfU - 1) ( 37 ) 

The optimality condition from dH/dlb = 0 yields 


1 


b 


T(r) 


Pv 


(38) 


The sufficient condition for the optimality condition (14) in this case is cPH/dll - 2pl 3 < 0 
with I 3 being an 3 x 3 identity matrix. We have p < 0, hence the well-known result that the 
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optimal thrust direction must be aligned with that of pv The costate equations (12) and 
(13) now become 


Pr = Pv 
PV = “Pr 


(39) 

(40) 


The costate equations have closed-form solution of 

cos r/3 sin rlz 


Pv(t) 


. -PrCO . 



Pv 0 

= n(r) 

PVo 

Pm 


P r o . 


where pvb and p ro are the (unknown) initial conditions for the costate. Define 

I c (r) = f l pv (C)cosCT(C)dC = f i c (CMC 
Jo Jo 

i s t) = f\ v (o^cnodc= r U<w 

Jo Jo 

It can be easily verified that the state equations have the solution of 


r(r) 

V(r) 


= U T ) 


ro 

V 0 


+ r(r) 


Ur) 

I.W 


(41) 

(42) 

(43) 

(44) 


where 


r (r) 


sin r/3 —cos r/3 
cos r/3 sin r/3 


(45) 


The integrals I c and I 3 can be evaluated by a numerical quadrature scheme. Calise et al 
use the Simpson’s rule in Ref. 7 We decide to use the Milne’s rule because it only increases 
computation by a small margin yet offers a combined four-orders of magnitude higher accu- 
racy for a 200-second burn. Let 6 = T togo /4 where Ttogo is the dimensionless time-to-go till 
MECO, yet to be determined. With Milne’s rule, we have 


I» (7*090) = hm + 321,(5) + 12ii(25) -1- 32ii(35) + 7i<(45)] , t = c, s (46) 


Note that we have left T vac inside the integrals (as a part of T(r)) because T vac will be 
time-varying (throttled down) when a thrust acceleration limit is imposed and becomes 
active. With the thnist integrals given in (46) and the costate given in (41) as function 
of its initial conditions, the state is found in closed-form from Eq. (45). Therefore the 
final state and costate are explicitly functions of pv„> Pm and Ttogo • Consequently, the total 6 
terminal conditions (10) and (18) are functions of the 7 unknowns pv„, Pr 0 and Ttogo ■ The 7th 
condition is from condition (17). For minimum-time problem, the costate can be scaled by 
arbitrary positive constant without changing any necessary conditions for the optimal control 
problem. And it will be shown that for Keplerian orbit insertion as the terminal conditions 
in (10), the condition (17) is automatically satisfied. Therefore the problem can actually be 
reduced to a six-unknown problem by requiring, for instance, that ||p 0 || = || (Pvj, Pr 0 )H = 
One of the six components of po is determined by the other five. But this still leaves us the 
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ambiguity of determining the sign of that component of po- We opt to avoid this problem 
by still treating the problem as a seven-unknown problem, and adding a trivial condition 

HpC^/)II = 1 

From the costate equations (39) and (40) , it can be shown by simple differentiation that 

d llP( T )H = o (48) 

dr 

Therefore ||p(r)|| =constant. The condition (47) thus will always be trivially satisfied if we 
scale po to have ||p 0 || = 1 (and we must, in order for the condition (47) to be met). 

For minimum-time problem where J = (f> = t/, the condition (17) be comes H\t { — 1 * 
This condition is equivalent to H\ tj > 0 because the costate p can always be scaled by 
a positive constant to achieve H\ tf = 1 if H \ tj > 0 before scaling. But this condition is 
automatically satisfied for Keplerian orbit insertion if other necessary conditions are satisfied. 
For Keplerian orbit insertion, the terminal conditions (10) represent the conditions on orbital 
elements which are constant if the engine is cut off at this point. Therefore 

<^( x /) = dfr( x /) ( v / ^ = o (49) 

dr dx f \~ T fJ 

where — r/ = g(r/) by the linear gravity approximation. Premultiply the above equation 
with the constant multiplier vector v in Eqs. (15-16) and use these transversality conditions. 

- T ^(^)=P^ V /-Pv, r / = 0 < 5 °) 

Use the expression of the Hamiltonian (37), l t = l Pv , and the above relationship, we arrive 
at 

H\t f =T(T f )\\p Vf \\>0 


In summary, the optimal vacuum ascent guidance problem becomes a root-finding prob- 
lem with seven unknowns (po and Tt ogo )> six constraints (10) and (18) plus one easy 
constraint (47). Through the use of quadrature (46), all the final state x/ and costate p / are 
explicit functions of the seven unknowns. The modified Newton method again works very 
well and the convergence occurs rapidly with almost any initial guesses that do not result in 
totally wrong initial thrust direction. 
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